function [dc, rc] = cfdFactorizeILU(ac,anb,bc,cconn)
%--------------------------------------------------------------------------
%
%  Written by the CFD Group @ AUB, Fall 2018
%  Contact us at: cfd@aub.edu.lb
%==========================================================================
% Routine Description:
%   Incomplete Lower Upper Factorization
%--------------------------------------------------------------------------

numberOfElements = length(ac);

dc = zeros(size(ac));
rc = zeros(size(ac));

for i1=1:numberOfElements
    dc(i1) = ac(i1); % dc就是对角线列
end

% 本质上在遍历所有的element，对于每个element i1（全局索引）
% 遍历i1的所有邻居element nb（全局索引jj1，在i1的邻居表内的局部索引为j1_）
% 在nb的邻居表内通过检查全局索引找到与i1对应的局部索引k1_
% 之所以需要局部索引是因为anb是将稀疏数据紧致存储的
for i1=1:numberOfElements
    dc(i1) = 1.0/dc(i1);
    rc(i1) = bc(i1);
    
    i1NbList = cconn{i1};
    i1NNb = length(i1NbList);
    
    if(i1~=numberOfElements-1)
        % loop over neighbours of iElement
        j1_ = 1;% element i1的邻居单元的局部索引
        while(j1_<=i1NNb)
            jj1 = i1NbList(j1_); % element i1的邻居element的全局索引为jj1
            % for all neighbour j > i do
            if((jj1>i1) && (jj1<=numberOfElements))
                j1NbList = cconn{jj1};% 邻居element jj1的邻居索引列表j1NbList
                j1NNb = length(j1NbList);
                i1_= 0;
                k1 = -1;% element i1的邻居单元（全局jj1，相对i1的局部j1_）的邻居单元的局部索引
                % find _i index to get A[j][_i] 找到k1 == i1为止，即从邻居
                while((i1_<=j1NNb) && (k1 ~= i1))
                    i1_ = i1_ + 1;
                    k1 = j1NbList(i1_);
                end
                % Compute A[j][i]*D[i]*A[i][j] 注意第一个下标用的是局部索引
                if(k1 == i1) % 执行到这里时dc(i1)已经执行了倒数，而dc(jj1)则还没有
                    dc(jj1) = dc(jj1) - anb{jj1}(i1_)*dc(i1)*anb{i1}(j1_);
                else
                    disp('the index for i in j is not found');
                end
            end
            j1_ = j1_ + 1;
        end
    end
end